<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN">
<html>
 <head>
  <title>slarre.f</title>
 <meta name="generator" content="emacs 21.3.1; htmlfontify 0.20">
<style type="text/css"><!-- 
body { background: rgb(255, 255, 255);  color: rgb(0, 0, 0);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: none; }
span.default   { background: rgb(255, 255, 255);  color: rgb(0, 0, 0);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: none; }
span.default a { background: rgb(255, 255, 255);  color: rgb(0, 0, 0);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: underline; }
span.string   { color: rgb(188, 143, 143);  background: rgb(255, 255, 255);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: none; }
span.string a { color: rgb(188, 143, 143);  background: rgb(255, 255, 255);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: underline; }
span.comment   { color: rgb(178, 34, 34);  background: rgb(255, 255, 255);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: none; }
span.comment a { color: rgb(178, 34, 34);  background: rgb(255, 255, 255);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: underline; }
 --></style>

 </head>
  <body>

<pre>
      SUBROUTINE <a name="SLARRE.1"></a><a href="slarre.f.html#SLARRE.1">SLARRE</a>( RANGE, N, VL, VU, IL, IU, D, E, E2,
     $                    RTOL1, RTOL2, SPLTOL, NSPLIT, ISPLIT, M,
     $                    W, WERR, WGAP, IBLOCK, INDEXW, GERS, PIVMIN,
     $                    WORK, IWORK, INFO )
      IMPLICIT NONE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  -- LAPACK auxiliary routine (version 3.1) --
</span><span class="comment">*</span><span class="comment">     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
</span><span class="comment">*</span><span class="comment">     November 2006
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     .. Scalar Arguments ..
</span>      CHARACTER          RANGE
      INTEGER            IL, INFO, IU, M, N, NSPLIT
      REAL              PIVMIN, RTOL1, RTOL2, SPLTOL, VL, VU
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Array Arguments ..
</span>      INTEGER            IBLOCK( * ), ISPLIT( * ), IWORK( * ),
     $                   INDEXW( * )
      REAL               D( * ), E( * ), E2( * ), GERS( * ),
     $                   W( * ),WERR( * ), WGAP( * ), WORK( * )
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  Purpose
</span><span class="comment">*</span><span class="comment">  =======
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  To find the desired eigenvalues of a given real symmetric
</span><span class="comment">*</span><span class="comment">  tridiagonal matrix T, <a name="SLARRE.27"></a><a href="slarre.f.html#SLARRE.1">SLARRE</a> sets any &quot;small&quot; off-diagonal
</span><span class="comment">*</span><span class="comment">  elements to zero, and for each unreduced block T_i, it finds
</span><span class="comment">*</span><span class="comment">  (a) a suitable shift at one end of the block's spectrum,
</span><span class="comment">*</span><span class="comment">  (b) the base representation, T_i - sigma_i I = L_i D_i L_i^T, and
</span><span class="comment">*</span><span class="comment">  (c) eigenvalues of each L_i D_i L_i^T.
</span><span class="comment">*</span><span class="comment">  The representations and eigenvalues found are then used by
</span><span class="comment">*</span><span class="comment">  <a name="SSTEMR.33"></a><a href="sstemr.f.html#SSTEMR.1">SSTEMR</a> to compute the eigenvectors of T.
</span><span class="comment">*</span><span class="comment">  The accuracy varies depending on whether bisection is used to
</span><span class="comment">*</span><span class="comment">  find a few eigenvalues or the dqds algorithm (subroutine <a name="SLASQ2.35"></a><a href="slasq2.f.html#SLASQ2.1">SLASQ2</a>) to
</span><span class="comment">*</span><span class="comment">  conpute all and then discard any unwanted one.
</span><span class="comment">*</span><span class="comment">  As an added benefit, <a name="SLARRE.37"></a><a href="slarre.f.html#SLARRE.1">SLARRE</a> also outputs the n
</span><span class="comment">*</span><span class="comment">  Gerschgorin intervals for the matrices L_i D_i L_i^T.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  Arguments
</span><span class="comment">*</span><span class="comment">  =========
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  RANGE   (input) CHARACTER
</span><span class="comment">*</span><span class="comment">          = 'A': (&quot;All&quot;)   all eigenvalues will be found.
</span><span class="comment">*</span><span class="comment">          = 'V': (&quot;Value&quot;) all eigenvalues in the half-open interval
</span><span class="comment">*</span><span class="comment">                           (VL, VU] will be found.
</span><span class="comment">*</span><span class="comment">          = 'I': (&quot;Index&quot;) the IL-th through IU-th eigenvalues (of the
</span><span class="comment">*</span><span class="comment">                           entire matrix) will be found.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  N       (input) INTEGER
</span><span class="comment">*</span><span class="comment">          The order of the matrix. N &gt; 0.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  VL      (input/output) REAL            
</span><span class="comment">*</span><span class="comment">  VU      (input/output) REAL            
</span><span class="comment">*</span><span class="comment">          If RANGE='V', the lower and upper bounds for the eigenvalues.
</span><span class="comment">*</span><span class="comment">          Eigenvalues less than or equal to VL, or greater than VU,
</span><span class="comment">*</span><span class="comment">          will not be returned.  VL &lt; VU.
</span><span class="comment">*</span><span class="comment">          If RANGE='I' or ='A', <a name="SLARRE.58"></a><a href="slarre.f.html#SLARRE.1">SLARRE</a> computes bounds on the desired
</span><span class="comment">*</span><span class="comment">          part of the spectrum.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  IL      (input) INTEGER
</span><span class="comment">*</span><span class="comment">  IU      (input) INTEGER
</span><span class="comment">*</span><span class="comment">          If RANGE='I', the indices (in ascending order) of the
</span><span class="comment">*</span><span class="comment">          smallest and largest eigenvalues to be returned.
</span><span class="comment">*</span><span class="comment">          1 &lt;= IL &lt;= IU &lt;= N.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  D       (input/output) REAL             array, dimension (N)
</span><span class="comment">*</span><span class="comment">          On entry, the N diagonal elements of the tridiagonal
</span><span class="comment">*</span><span class="comment">          matrix T.
</span><span class="comment">*</span><span class="comment">          On exit, the N diagonal elements of the diagonal
</span><span class="comment">*</span><span class="comment">          matrices D_i.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  E       (input/output) REAL             array, dimension (N)
</span><span class="comment">*</span><span class="comment">          On entry, the first (N-1) entries contain the subdiagonal
</span><span class="comment">*</span><span class="comment">          elements of the tridiagonal matrix T; E(N) need not be set.
</span><span class="comment">*</span><span class="comment">          On exit, E contains the subdiagonal elements of the unit
</span><span class="comment">*</span><span class="comment">          bidiagonal matrices L_i. The entries E( ISPLIT( I ) ),
</span><span class="comment">*</span><span class="comment">          1 &lt;= I &lt;= NSPLIT, contain the base points sigma_i on output.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  E2      (input/output) REAL             array, dimension (N)
</span><span class="comment">*</span><span class="comment">          On entry, the first (N-1) entries contain the SQUARES of the
</span><span class="comment">*</span><span class="comment">          subdiagonal elements of the tridiagonal matrix T;
</span><span class="comment">*</span><span class="comment">          E2(N) need not be set.
</span><span class="comment">*</span><span class="comment">          On exit, the entries E2( ISPLIT( I ) ),
</span><span class="comment">*</span><span class="comment">          1 &lt;= I &lt;= NSPLIT, have been set to zero
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  RTOL1   (input) REAL            
</span><span class="comment">*</span><span class="comment">  RTOL2   (input) REAL            
</span><span class="comment">*</span><span class="comment">           Parameters for bisection.
</span><span class="comment">*</span><span class="comment">           An interval [LEFT,RIGHT] has converged if
</span><span class="comment">*</span><span class="comment">           RIGHT-LEFT.LT.MAX( RTOL1*GAP, RTOL2*MAX(|LEFT|,|RIGHT|) )
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  SPLTOL (input) REAL            
</span><span class="comment">*</span><span class="comment">          The threshold for splitting.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  NSPLIT  (output) INTEGER
</span><span class="comment">*</span><span class="comment">          The number of blocks T splits into. 1 &lt;= NSPLIT &lt;= N.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  ISPLIT  (output) INTEGER array, dimension (N)
</span><span class="comment">*</span><span class="comment">          The splitting points, at which T breaks up into blocks.
</span><span class="comment">*</span><span class="comment">          The first block consists of rows/columns 1 to ISPLIT(1),
</span><span class="comment">*</span><span class="comment">          the second of rows/columns ISPLIT(1)+1 through ISPLIT(2),
</span><span class="comment">*</span><span class="comment">          etc., and the NSPLIT-th consists of rows/columns
</span><span class="comment">*</span><span class="comment">          ISPLIT(NSPLIT-1)+1 through ISPLIT(NSPLIT)=N.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  M       (output) INTEGER
</span><span class="comment">*</span><span class="comment">          The total number of eigenvalues (of all L_i D_i L_i^T)
</span><span class="comment">*</span><span class="comment">          found.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  W       (output) REAL             array, dimension (N)
</span><span class="comment">*</span><span class="comment">          The first M elements contain the eigenvalues. The
</span><span class="comment">*</span><span class="comment">          eigenvalues of each of the blocks, L_i D_i L_i^T, are
</span><span class="comment">*</span><span class="comment">          sorted in ascending order ( <a name="SLARRE.113"></a><a href="slarre.f.html#SLARRE.1">SLARRE</a> may use the
</span><span class="comment">*</span><span class="comment">          remaining N-M elements as workspace).
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  WERR    (output) REAL             array, dimension (N)
</span><span class="comment">*</span><span class="comment">          The error bound on the corresponding eigenvalue in W.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  WGAP    (output) REAL             array, dimension (N)
</span><span class="comment">*</span><span class="comment">          The separation from the right neighbor eigenvalue in W.
</span><span class="comment">*</span><span class="comment">          The gap is only with respect to the eigenvalues of the same block
</span><span class="comment">*</span><span class="comment">          as each block has its own representation tree.
</span><span class="comment">*</span><span class="comment">          Exception: at the right end of a block we store the left gap
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  IBLOCK  (output) INTEGER array, dimension (N)
</span><span class="comment">*</span><span class="comment">          The indices of the blocks (submatrices) associated with the
</span><span class="comment">*</span><span class="comment">          corresponding eigenvalues in W; IBLOCK(i)=1 if eigenvalue
</span><span class="comment">*</span><span class="comment">          W(i) belongs to the first block from the top, =2 if W(i)
</span><span class="comment">*</span><span class="comment">          belongs to the second block, etc.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  INDEXW  (output) INTEGER array, dimension (N)
</span><span class="comment">*</span><span class="comment">          The indices of the eigenvalues within each block (submatrix);
</span><span class="comment">*</span><span class="comment">          for example, INDEXW(i)= 10 and IBLOCK(i)=2 imply that the
</span><span class="comment">*</span><span class="comment">          i-th eigenvalue W(i) is the 10-th eigenvalue in block 2
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  GERS    (output) REAL             array, dimension (2*N)
</span><span class="comment">*</span><span class="comment">          The N Gerschgorin intervals (the i-th Gerschgorin interval
</span><span class="comment">*</span><span class="comment">          is (GERS(2*i-1), GERS(2*i)).
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  PIVMIN  (output) DOUBLE PRECISION
</span><span class="comment">*</span><span class="comment">          The minimum pivot in the Sturm sequence for T.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  WORK    (workspace) REAL             array, dimension (6*N)
</span><span class="comment">*</span><span class="comment">          Workspace.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  IWORK   (workspace) INTEGER array, dimension (5*N)
</span><span class="comment">*</span><span class="comment">          Workspace.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  INFO    (output) INTEGER
</span><span class="comment">*</span><span class="comment">          = 0:  successful exit
</span><span class="comment">*</span><span class="comment">          &gt; 0:  A problem occured in <a name="SLARRE.151"></a><a href="slarre.f.html#SLARRE.1">SLARRE</a>.
</span><span class="comment">*</span><span class="comment">          &lt; 0:  One of the called subroutines signaled an internal problem.
</span><span class="comment">*</span><span class="comment">                Needs inspection of the corresponding parameter IINFO
</span><span class="comment">*</span><span class="comment">                for further information.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">          =-1:  Problem in <a name="SLARRD.156"></a><a href="slarrd.f.html#SLARRD.1">SLARRD</a>.
</span><span class="comment">*</span><span class="comment">          = 2:  No base representation could be found in MAXTRY iterations.
</span><span class="comment">*</span><span class="comment">                Increasing MAXTRY and recompilation might be a remedy.
</span><span class="comment">*</span><span class="comment">          =-3:  Problem in <a name="SLARRB.159"></a><a href="slarrb.f.html#SLARRB.1">SLARRB</a> when computing the refined root
</span><span class="comment">*</span><span class="comment">                representation for <a name="SLASQ2.160"></a><a href="slasq2.f.html#SLASQ2.1">SLASQ2</a>.
</span><span class="comment">*</span><span class="comment">          =-4:  Problem in <a name="SLARRB.161"></a><a href="slarrb.f.html#SLARRB.1">SLARRB</a> when preforming bisection on the
</span><span class="comment">*</span><span class="comment">                desired part of the spectrum.
</span><span class="comment">*</span><span class="comment">          =-5:  Problem in <a name="SLASQ2.163"></a><a href="slasq2.f.html#SLASQ2.1">SLASQ2</a>.
</span><span class="comment">*</span><span class="comment">          =-6:  Problem in <a name="SLASQ2.164"></a><a href="slasq2.f.html#SLASQ2.1">SLASQ2</a>.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  Further Details
</span><span class="comment">*</span><span class="comment">  The base representations are required to suffer very little
</span><span class="comment">*</span><span class="comment">  element growth and consequently define all their eigenvalues to
</span><span class="comment">*</span><span class="comment">  high relative accuracy.
</span><span class="comment">*</span><span class="comment">  ===============
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  Based on contributions by
</span><span class="comment">*</span><span class="comment">     Beresford Parlett, University of California, Berkeley, USA
</span><span class="comment">*</span><span class="comment">     Jim Demmel, University of California, Berkeley, USA
</span><span class="comment">*</span><span class="comment">     Inderjit Dhillon, University of Texas, Austin, USA
</span><span class="comment">*</span><span class="comment">     Osni Marques, LBNL/NERSC, USA
</span><span class="comment">*</span><span class="comment">     Christof Voemel, University of California, Berkeley, USA
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  =====================================================================
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     .. Parameters ..
</span>      REAL               FAC, FOUR, FOURTH, FUDGE, HALF, HNDRD,
     $                   MAXGROWTH, ONE, PERT, TWO, ZERO
      PARAMETER          ( ZERO = 0.0E0, ONE = 1.0E0,
     $                     TWO = 2.0E0, FOUR=4.0E0,
     $                     HNDRD = 100.0E0,
     $                     PERT = 4.0E0,
     $                     HALF = ONE/TWO, FOURTH = ONE/FOUR, FAC= HALF,
     $                     MAXGROWTH = 64.0E0, FUDGE = 2.0E0 )
      INTEGER            MAXTRY, ALLRNG, INDRNG, VALRNG
      PARAMETER          ( MAXTRY = 6, ALLRNG = 1, INDRNG = 2,
     $                     VALRNG = 3 )
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Local Scalars ..
</span>      LOGICAL            FORCEB, NOREP, USEDQD
      INTEGER            CNT, CNT1, CNT2, I, IBEGIN, IDUM, IEND, IINFO,
     $                   IN, INDL, INDU, IRANGE, J, JBLK, MB, MM,
     $                   WBEGIN, WEND
      REAL               AVGAP, BSRTOL, CLWDTH, DMAX, DPIVOT, EABS,
     $                   EMAX, EOLD, EPS, GL, GU, ISLEFT, ISRGHT, RTL,
     $                   RTOL, S1, S2, SAFMIN, SGNDEF, SIGMA, SPDIAM,
     $                   TAU, TMP, TMP1


<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Local Arrays ..
</span>      INTEGER            ISEED( 4 )
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. External Functions ..
</span>      LOGICAL            <a name="LSAME.210"></a><a href="lsame.f.html#LSAME.1">LSAME</a>
      REAL                        <a name="SLAMCH.211"></a><a href="slamch.f.html#SLAMCH.1">SLAMCH</a>
      EXTERNAL           <a name="SLAMCH.212"></a><a href="slamch.f.html#SLAMCH.1">SLAMCH</a>, <a name="LSAME.212"></a><a href="lsame.f.html#LSAME.1">LSAME</a>

<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. External Subroutines ..
</span>      EXTERNAL           SCOPY, <a name="SLARNV.216"></a><a href="slarnv.f.html#SLARNV.1">SLARNV</a>, <a name="SLARRA.216"></a><a href="slarra.f.html#SLARRA.1">SLARRA</a>, <a name="SLARRB.216"></a><a href="slarrb.f.html#SLARRB.1">SLARRB</a>, <a name="SLARRC.216"></a><a href="slarrc.f.html#SLARRC.1">SLARRC</a>, <a name="SLARRD.216"></a><a href="slarrd.f.html#SLARRD.1">SLARRD</a>,
     $                   <a name="SLASQ2.217"></a><a href="slasq2.f.html#SLASQ2.1">SLASQ2</a>
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Intrinsic Functions ..
</span>      INTRINSIC          ABS, MAX, MIN

<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Executable Statements ..
</span><span class="comment">*</span><span class="comment">
</span>
      INFO = 0

<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Decode RANGE
</span><span class="comment">*</span><span class="comment">
</span>      IF( <a name="LSAME.231"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( RANGE, <span class="string">'A'</span> ) ) THEN
         IRANGE = ALLRNG
      ELSE IF( <a name="LSAME.233"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( RANGE, <span class="string">'V'</span> ) ) THEN
         IRANGE = VALRNG
      ELSE IF( <a name="LSAME.235"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( RANGE, <span class="string">'I'</span> ) ) THEN
         IRANGE = INDRNG
      END IF

      M = 0

<span class="comment">*</span><span class="comment">     Get machine constants
</span>      SAFMIN = <a name="SLAMCH.242"></a><a href="slamch.f.html#SLAMCH.1">SLAMCH</a>( <span class="string">'S'</span> )
      EPS = <a name="SLAMCH.243"></a><a href="slamch.f.html#SLAMCH.1">SLAMCH</a>( <span class="string">'P'</span> )

<span class="comment">*</span><span class="comment">     Set parameters
</span>      RTL = HNDRD*EPS
<span class="comment">*</span><span class="comment">     If one were ever to ask for less initial precision in BSRTOL,
</span><span class="comment">*</span><span class="comment">     one should keep in mind that for the subset case, the extremal
</span><span class="comment">*</span><span class="comment">     eigenvalues must be at least as accurate as the current setting
</span><span class="comment">*</span><span class="comment">     (eigenvalues in the middle need not as much accuracy)
</span>      BSRTOL = SQRT(EPS)*(0.5E-3)

<span class="comment">*</span><span class="comment">     Treat case of 1x1 matrix for quick return
</span>      IF( N.EQ.1 ) THEN
         IF( (IRANGE.EQ.ALLRNG).OR.
     $       ((IRANGE.EQ.VALRNG).AND.(D(1).GT.VL).AND.(D(1).LE.VU)).OR.
     $       ((IRANGE.EQ.INDRNG).AND.(IL.EQ.1).AND.(IU.EQ.1)) ) THEN
            M = 1
            W(1) = D(1)
<span class="comment">*</span><span class="comment">           The computation error of the eigenvalue is zero
</span>            WERR(1) = ZERO
            WGAP(1) = ZERO
            IBLOCK( 1 ) = 1
            INDEXW( 1 ) = 1
            GERS(1) = D( 1 )
            GERS(2) = D( 1 )
         ENDIF
<span class="comment">*</span><span class="comment">        store the shift for the initial RRR, which is zero in this case
</span>         E(1) = ZERO
         RETURN
      END IF

<span class="comment">*</span><span class="comment">     General case: tridiagonal matrix of order &gt; 1
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Init WERR, WGAP. Compute Gerschgorin intervals and spectral diameter.
</span><span class="comment">*</span><span class="comment">     Compute maximum off-diagonal entry and pivmin.
</span>      GL = D(1)
      GU = D(1)
      EOLD = ZERO
      EMAX = ZERO
      E(N) = ZERO
      DO 5 I = 1,N
         WERR(I) = ZERO
         WGAP(I) = ZERO
         EABS = ABS( E(I) )
         IF( EABS .GE. EMAX ) THEN
            EMAX = EABS
         END IF
         TMP1 = EABS + EOLD
         GERS( 2*I-1) = D(I) - TMP1
         GL =  MIN( GL, GERS( 2*I - 1))
         GERS( 2*I ) = D(I) + TMP1
         GU = MAX( GU, GERS(2*I) )
         EOLD  = EABS
 5    CONTINUE
<span class="comment">*</span><span class="comment">     The minimum pivot allowed in the Sturm sequence for T
</span>      PIVMIN = SAFMIN * MAX( ONE, EMAX**2 )
<span class="comment">*</span><span class="comment">     Compute spectral diameter. The Gerschgorin bounds give an
</span><span class="comment">*</span><span class="comment">     estimate that is wrong by at most a factor of SQRT(2)
</span>      SPDIAM = GU - GL

<span class="comment">*</span><span class="comment">     Compute splitting points
</span>      CALL <a name="SLARRA.303"></a><a href="slarra.f.html#SLARRA.1">SLARRA</a>( N, D, E, E2, SPLTOL, SPDIAM,
     $                    NSPLIT, ISPLIT, IINFO )

<span class="comment">*</span><span class="comment">     Can force use of bisection instead of faster DQDS.
</span><span class="comment">*</span><span class="comment">     Option left in the code for future multisection work.
</span>      FORCEB = .FALSE.

      IF( (IRANGE.EQ.ALLRNG) .AND. (.NOT. FORCEB) ) THEN
<span class="comment">*</span><span class="comment">        Set interval [VL,VU] that contains all eigenvalues
</span>         VL = GL
         VU = GU
      ELSE
<span class="comment">*</span><span class="comment">        We call <a name="SLARRD.315"></a><a href="slarrd.f.html#SLARRD.1">SLARRD</a> to find crude approximations to the eigenvalues
</span><span class="comment">*</span><span class="comment">        in the desired range. In case IRANGE = INDRNG, we also obtain the
</span><span class="comment">*</span><span class="comment">        interval (VL,VU] that contains all the wanted eigenvalues.
</span><span class="comment">*</span><span class="comment">        An interval [LEFT,RIGHT] has converged if
</span><span class="comment">*</span><span class="comment">        RIGHT-LEFT.LT.RTOL*MAX(ABS(LEFT),ABS(RIGHT))
</span><span class="comment">*</span><span class="comment">        <a name="SLARRD.320"></a><a href="slarrd.f.html#SLARRD.1">SLARRD</a> needs a WORK of size 4*N, IWORK of size 3*N
</span>         CALL <a name="SLARRD.321"></a><a href="slarrd.f.html#SLARRD.1">SLARRD</a>( RANGE, <span class="string">'B'</span>, N, VL, VU, IL, IU, GERS,
     $                    BSRTOL, D, E, E2, PIVMIN, NSPLIT, ISPLIT,
     $                    MM, W, WERR, VL, VU, IBLOCK, INDEXW,
     $                    WORK, IWORK, IINFO )
         IF( IINFO.NE.0 ) THEN
            INFO = -1
            RETURN
         ENDIF
<span class="comment">*</span><span class="comment">        Make sure that the entries M+1 to N in W, WERR, IBLOCK, INDEXW are 0
</span>         DO 14 I = MM+1,N
            W( I ) = ZERO
            WERR( I ) = ZERO
            IBLOCK( I ) = 0
            INDEXW( I ) = 0
 14      CONTINUE
      END IF


<span class="comment">*</span><span class="comment">**
</span><span class="comment">*</span><span class="comment">     Loop over unreduced blocks
</span>      IBEGIN = 1
      WBEGIN = 1
      DO 170 JBLK = 1, NSPLIT
         IEND = ISPLIT( JBLK )
         IN = IEND - IBEGIN + 1

<span class="comment">*</span><span class="comment">        1 X 1 block
</span>         IF( IN.EQ.1 ) THEN
            IF( (IRANGE.EQ.ALLRNG).OR.( (IRANGE.EQ.VALRNG).AND.
     $         ( D( IBEGIN ).GT.VL ).AND.( D( IBEGIN ).LE.VU ) )
     $        .OR. ( (IRANGE.EQ.INDRNG).AND.(IBLOCK(WBEGIN).EQ.JBLK))
     $        ) THEN
               M = M + 1
               W( M ) = D( IBEGIN )
               WERR(M) = ZERO
<span class="comment">*</span><span class="comment">              The gap for a single block doesn't matter for the later
</span><span class="comment">*</span><span class="comment">              algorithm and is assigned an arbitrary large value
</span>               WGAP(M) = ZERO
               IBLOCK( M ) = JBLK
               INDEXW( M ) = 1
               WBEGIN = WBEGIN + 1
            ENDIF
<span class="comment">*</span><span class="comment">           E( IEND ) holds the shift for the initial RRR
</span>            E( IEND ) = ZERO
            IBEGIN = IEND + 1
            GO TO 170
         END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Blocks of size larger than 1x1
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        E( IEND ) will hold the shift for the initial RRR, for now set it =0
</span>         E( IEND ) = ZERO
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Find local outer bounds GL,GU for the block
</span>         GL = D(IBEGIN)
         GU = D(IBEGIN)
         DO 15 I = IBEGIN , IEND
            GL = MIN( GERS( 2*I-1 ), GL )
            GU = MAX( GERS( 2*I ), GU )
 15      CONTINUE
         SPDIAM = GU - GL

         IF(.NOT. ((IRANGE.EQ.ALLRNG).AND.(.NOT.FORCEB)) ) THEN
<span class="comment">*</span><span class="comment">           Count the number of eigenvalues in the current block.
</span>            MB = 0
            DO 20 I = WBEGIN,MM
               IF( IBLOCK(I).EQ.JBLK ) THEN
                  MB = MB+1
               ELSE
                  GOTO 21
               ENDIF
 20         CONTINUE
 21         CONTINUE

            IF( MB.EQ.0) THEN
<span class="comment">*</span><span class="comment">              No eigenvalue in the current block lies in the desired range
</span><span class="comment">*</span><span class="comment">              E( IEND ) holds the shift for the initial RRR
</span>               E( IEND ) = ZERO
               IBEGIN = IEND + 1
               GO TO 170
            ELSE

<span class="comment">*</span><span class="comment">              Decide whether dqds or bisection is more efficient
</span>               USEDQD = ( (MB .GT. FAC*IN) .AND. (.NOT.FORCEB) )
               WEND = WBEGIN + MB - 1
<span class="comment">*</span><span class="comment">              Calculate gaps for the current block
</span><span class="comment">*</span><span class="comment">              In later stages, when representations for individual
</span><span class="comment">*</span><span class="comment">              eigenvalues are different, we use SIGMA = E( IEND ).
</span>               SIGMA = ZERO
               DO 30 I = WBEGIN, WEND - 1
                  WGAP( I ) = MAX( ZERO,
     $                        W(I+1)-WERR(I+1) - (W(I)+WERR(I)) )
 30            CONTINUE
               WGAP( WEND ) = MAX( ZERO,
     $                     VU - SIGMA - (W( WEND )+WERR( WEND )))
<span class="comment">*</span><span class="comment">              Find local index of the first and last desired evalue.
</span>               INDL = INDEXW(WBEGIN)
               INDU = INDEXW( WEND )
            ENDIF
         ENDIF
         IF(( (IRANGE.EQ.ALLRNG) .AND. (.NOT. FORCEB) ).OR.USEDQD) THEN
<span class="comment">*</span><span class="comment">           Case of DQDS
</span><span class="comment">*</span><span class="comment">           Find approximations to the extremal eigenvalues of the block
</span>            CALL <a name="SLARRK.424"></a><a href="slarrk.f.html#SLARRK.1">SLARRK</a>( IN, 1, GL, GU, D(IBEGIN),
     $               E2(IBEGIN), PIVMIN, RTL, TMP, TMP1, IINFO )
            IF( IINFO.NE.0 ) THEN
               INFO = -1
               RETURN
            ENDIF
            ISLEFT = MAX(GL, TMP - TMP1
     $               - HNDRD * EPS* ABS(TMP - TMP1))

            CALL <a name="SLARRK.433"></a><a href="slarrk.f.html#SLARRK.1">SLARRK</a>( IN, IN, GL, GU, D(IBEGIN),
     $               E2(IBEGIN), PIVMIN, RTL, TMP, TMP1, IINFO )
            IF( IINFO.NE.0 ) THEN
               INFO = -1
               RETURN
            ENDIF
            ISRGHT = MIN(GU, TMP + TMP1
     $                 + HNDRD * EPS * ABS(TMP + TMP1))
<span class="comment">*</span><span class="comment">           Improve the estimate of the spectral diameter
</span>            SPDIAM = ISRGHT - ISLEFT
         ELSE
<span class="comment">*</span><span class="comment">           Case of bisection
</span><span class="comment">*</span><span class="comment">           Find approximations to the wanted extremal eigenvalues
</span>            ISLEFT = MAX(GL, W(WBEGIN) - WERR(WBEGIN)
     $                  - HNDRD * EPS*ABS(W(WBEGIN)- WERR(WBEGIN) ))
            ISRGHT = MIN(GU,W(WEND) + WERR(WEND)
     $                  + HNDRD * EPS * ABS(W(WEND)+ WERR(WEND)))
         ENDIF


<span class="comment">*</span><span class="comment">        Decide whether the base representation for the current block
</span><span class="comment">*</span><span class="comment">        L_JBLK D_JBLK L_JBLK^T = T_JBLK - sigma_JBLK I
</span><span class="comment">*</span><span class="comment">        should be on the left or the right end of the current block.
</span><span class="comment">*</span><span class="comment">        The strategy is to shift to the end which is &quot;more populated&quot;
</span><span class="comment">*</span><span class="comment">        Furthermore, decide whether to use DQDS for the computation of
</span><span class="comment">*</span><span class="comment">        the eigenvalue approximations at the end of <a name="SLARRE.458"></a><a href="slarre.f.html#SLARRE.1">SLARRE</a> or bisection.
</span><span class="comment">*</span><span class="comment">        dqds is chosen if all eigenvalues are desired or the number of
</span><span class="comment">*</span><span class="comment">        eigenvalues to be computed is large compared to the blocksize.
</span>         IF( ( IRANGE.EQ.ALLRNG ) .AND. (.NOT.FORCEB) ) THEN
<span class="comment">*</span><span class="comment">           If all the eigenvalues have to be computed, we use dqd
</span>            USEDQD = .TRUE.
<span class="comment">*</span><span class="comment">           INDL is the local index of the first eigenvalue to compute
</span>            INDL = 1
            INDU = IN
<span class="comment">*</span><span class="comment">           MB =  number of eigenvalues to compute
</span>            MB = IN
            WEND = WBEGIN + MB - 1
<span class="comment">*</span><span class="comment">           Define 1/4 and 3/4 points of the spectrum
</span>            S1 = ISLEFT + FOURTH * SPDIAM
            S2 = ISRGHT - FOURTH * SPDIAM
         ELSE
<span class="comment">*</span><span class="comment">           <a name="SLARRD.474"></a><a href="slarrd.f.html#SLARRD.1">SLARRD</a> has computed IBLOCK and INDEXW for each eigenvalue
</span><span class="comment">*</span><span class="comment">           approximation.
</span><span class="comment">*</span><span class="comment">           choose sigma
</span>            IF( USEDQD ) THEN
               S1 = ISLEFT + FOURTH * SPDIAM
               S2 = ISRGHT - FOURTH * SPDIAM
            ELSE
               TMP = MIN(ISRGHT,VU) -  MAX(ISLEFT,VL)
               S1 =  MAX(ISLEFT,VL) + FOURTH * TMP
               S2 =  MIN(ISRGHT,VU) - FOURTH * TMP
            ENDIF
         ENDIF

<span class="comment">*</span><span class="comment">        Compute the negcount at the 1/4 and 3/4 points
</span>         IF(MB.GT.1) THEN
            CALL <a name="SLARRC.489"></a><a href="slarrc.f.html#SLARRC.1">SLARRC</a>( <span class="string">'T'</span>, IN, S1, S2, D(IBEGIN),
     $                    E(IBEGIN), PIVMIN, CNT, CNT1, CNT2, IINFO)
         ENDIF

         IF(MB.EQ.1) THEN
            SIGMA = GL
            SGNDEF = ONE
         ELSEIF( CNT1 - INDL .GE. INDU - CNT2 ) THEN
            IF( ( IRANGE.EQ.ALLRNG ) .AND. (.NOT.FORCEB) ) THEN
               SIGMA = MAX(ISLEFT,GL)
            ELSEIF( USEDQD ) THEN
<span class="comment">*</span><span class="comment">              use Gerschgorin bound as shift to get pos def matrix
</span><span class="comment">*</span><span class="comment">              for dqds
</span>               SIGMA = ISLEFT
            ELSE
<span class="comment">*</span><span class="comment">              use approximation of the first desired eigenvalue of the
</span><span class="comment">*</span><span class="comment">              block as shift
</span>               SIGMA = MAX(ISLEFT,VL)
            ENDIF
            SGNDEF = ONE
         ELSE
            IF( ( IRANGE.EQ.ALLRNG ) .AND. (.NOT.FORCEB) ) THEN
               SIGMA = MIN(ISRGHT,GU)
            ELSEIF( USEDQD ) THEN
<span class="comment">*</span><span class="comment">              use Gerschgorin bound as shift to get neg def matrix
</span><span class="comment">*</span><span class="comment">              for dqds
</span>               SIGMA = ISRGHT
            ELSE
<span class="comment">*</span><span class="comment">              use approximation of the first desired eigenvalue of the
</span><span class="comment">*</span><span class="comment">              block as shift
</span>               SIGMA = MIN(ISRGHT,VU)
            ENDIF
            SGNDEF = -ONE
         ENDIF


<span class="comment">*</span><span class="comment">        An initial SIGMA has been chosen that will be used for computing
</span><span class="comment">*</span><span class="comment">        T - SIGMA I = L D L^T
</span><span class="comment">*</span><span class="comment">        Define the increment TAU of the shift in case the initial shift
</span><span class="comment">*</span><span class="comment">        needs to be refined to obtain a factorization with not too much
</span><span class="comment">*</span><span class="comment">        element growth.
</span>         IF( USEDQD ) THEN
<span class="comment">*</span><span class="comment">           The initial SIGMA was to the outer end of the spectrum
</span><span class="comment">*</span><span class="comment">           the matrix is definite and we need not retreat.
</span>            TAU = SPDIAM*EPS*N + TWO*PIVMIN
         ELSE
            IF(MB.GT.1) THEN
               CLWDTH = W(WEND) + WERR(WEND) - W(WBEGIN) - WERR(WBEGIN)
               AVGAP = ABS(CLWDTH / REAL(WEND-WBEGIN))
               IF( SGNDEF.EQ.ONE ) THEN
                  TAU = HALF*MAX(WGAP(WBEGIN),AVGAP)
                  TAU = MAX(TAU,WERR(WBEGIN))
               ELSE
                  TAU = HALF*MAX(WGAP(WEND-1),AVGAP)
                  TAU = MAX(TAU,WERR(WEND))
               ENDIF
            ELSE
               TAU = WERR(WBEGIN)
            ENDIF
         ENDIF
<span class="comment">*</span><span class="comment">
</span>         DO 80 IDUM = 1, MAXTRY
<span class="comment">*</span><span class="comment">           Compute L D L^T factorization of tridiagonal matrix T - sigma I.
</span><span class="comment">*</span><span class="comment">           Store D in WORK(1:IN), L in WORK(IN+1:2*IN), and reciprocals of
</span><span class="comment">*</span><span class="comment">           pivots in WORK(2*IN+1:3*IN)
</span>            DPIVOT = D( IBEGIN ) - SIGMA
            WORK( 1 ) = DPIVOT
            DMAX = ABS( WORK(1) )
            J = IBEGIN
            DO 70 I = 1, IN - 1
               WORK( 2*IN+I ) = ONE / WORK( I )
               TMP = E( J )*WORK( 2*IN+I )
               WORK( IN+I ) = TMP
               DPIVOT = ( D( J+1 )-SIGMA ) - TMP*E( J )
               WORK( I+1 ) = DPIVOT
               DMAX = MAX( DMAX, ABS(DPIVOT) )
               J = J + 1
 70         CONTINUE
<span class="comment">*</span><span class="comment">           check for element growth
</span>            IF( DMAX .GT. MAXGROWTH*SPDIAM ) THEN
               NOREP = .TRUE.
            ELSE
               NOREP = .FALSE.
            ENDIF
            IF( USEDQD .AND. .NOT.NOREP ) THEN
<span class="comment">*</span><span class="comment">              Ensure the definiteness of the representation
</span><span class="comment">*</span><span class="comment">              All entries of D (of L D L^T) must have the same sign
</span>               DO 71 I = 1, IN
                  TMP = SGNDEF*WORK( I )
                  IF( TMP.LT.ZERO ) NOREP = .TRUE.
 71            CONTINUE
            ENDIF
            IF(NOREP) THEN
<span class="comment">*</span><span class="comment">              Note that in the case of IRANGE=ALLRNG, we use the Gerschgorin
</span><span class="comment">*</span><span class="comment">              shift which makes the matrix definite. So we should end up
</span><span class="comment">*</span><span class="comment">              here really only in the case of IRANGE = VALRNG or INDRNG.
</span>               IF( IDUM.EQ.MAXTRY-1 ) THEN
                  IF( SGNDEF.EQ.ONE ) THEN
<span class="comment">*</span><span class="comment">                    The fudged Gerschgorin shift should succeed
</span>                     SIGMA =
     $                    GL - FUDGE*SPDIAM*EPS*N - FUDGE*TWO*PIVMIN
                  ELSE
                     SIGMA =
     $                    GU + FUDGE*SPDIAM*EPS*N + FUDGE*TWO*PIVMIN
                  END IF
               ELSE
                  SIGMA = SIGMA - SGNDEF * TAU
                  TAU = TWO * TAU
               END IF
            ELSE
<span class="comment">*</span><span class="comment">              an initial RRR is found
</span>               GO TO 83
            END IF
 80      CONTINUE
<span class="comment">*</span><span class="comment">        if the program reaches this point, no base representation could be
</span><span class="comment">*</span><span class="comment">        found in MAXTRY iterations.
</span>         INFO = 2
         RETURN

 83      CONTINUE
<span class="comment">*</span><span class="comment">        At this point, we have found an initial base representation
</span><span class="comment">*</span><span class="comment">        T - SIGMA I = L D L^T with not too much element growth.
</span><span class="comment">*</span><span class="comment">        Store the shift.
</span>         E( IEND ) = SIGMA
<span class="comment">*</span><span class="comment">        Store D and L.
</span>         CALL SCOPY( IN, WORK, 1, D( IBEGIN ), 1 )
         CALL SCOPY( IN-1, WORK( IN+1 ), 1, E( IBEGIN ), 1 )


         IF(MB.GT.1 ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Perturb each entry of the base representation by a small
</span><span class="comment">*</span><span class="comment">           (but random) relative amount to overcome difficulties with
</span><span class="comment">*</span><span class="comment">           glued matrices.
</span><span class="comment">*</span><span class="comment">
</span>            DO 122 I = 1, 4
               ISEED( I ) = 1
 122        CONTINUE

            CALL <a name="SLARNV.628"></a><a href="slarnv.f.html#SLARNV.1">SLARNV</a>(2, ISEED, 2*IN-1, WORK(1))
            DO 125 I = 1,IN-1
               D(IBEGIN+I-1) = D(IBEGIN+I-1)*(ONE+EPS*PERT*WORK(I))
               E(IBEGIN+I-1) = E(IBEGIN+I-1)*(ONE+EPS*PERT*WORK(IN+I))
 125        CONTINUE
            D(IEND) = D(IEND)*(ONE+EPS*FOUR*WORK(IN))
<span class="comment">*</span><span class="comment">
</span>         ENDIF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Don't update the Gerschgorin intervals because keeping track
</span><span class="comment">*</span><span class="comment">        of the updates would be too much work in <a name="SLARRV.638"></a><a href="slarrv.f.html#SLARRV.1">SLARRV</a>.
</span><span class="comment">*</span><span class="comment">        We update W instead and use it to locate the proper Gerschgorin
</span><span class="comment">*</span><span class="comment">        intervals.
</span>
<span class="comment">*</span><span class="comment">        Compute the required eigenvalues of L D L' by bisection or dqds
</span>         IF ( .NOT.USEDQD ) THEN
<span class="comment">*</span><span class="comment">           If <a name="SLARRD.644"></a><a href="slarrd.f.html#SLARRD.1">SLARRD</a> has been used, shift the eigenvalue approximations
</span><span class="comment">*</span><span class="comment">           according to their representation. This is necessary for
</span><span class="comment">*</span><span class="comment">           a uniform <a name="SLARRV.646"></a><a href="slarrv.f.html#SLARRV.1">SLARRV</a> since dqds computes eigenvalues of the
</span><span class="comment">*</span><span class="comment">           shifted representation. In <a name="SLARRV.647"></a><a href="slarrv.f.html#SLARRV.1">SLARRV</a>, W will always hold the
</span><span class="comment">*</span><span class="comment">           UNshifted eigenvalue approximation.
</span>            DO 134 J=WBEGIN,WEND
               W(J) = W(J) - SIGMA
               WERR(J) = WERR(J) + ABS(W(J)) * EPS
 134        CONTINUE
<span class="comment">*</span><span class="comment">           call <a name="SLARRB.653"></a><a href="slarrb.f.html#SLARRB.1">SLARRB</a> to reduce eigenvalue error of the approximations
</span><span class="comment">*</span><span class="comment">           from <a name="SLARRD.654"></a><a href="slarrd.f.html#SLARRD.1">SLARRD</a>
</span>            DO 135 I = IBEGIN, IEND-1
               WORK( I ) = D( I ) * E( I )**2
 135        CONTINUE
<span class="comment">*</span><span class="comment">           use bisection to find EV from INDL to INDU
</span>            CALL <a name="SLARRB.659"></a><a href="slarrb.f.html#SLARRB.1">SLARRB</a>(IN, D(IBEGIN), WORK(IBEGIN),
     $                  INDL, INDU, RTOL1, RTOL2, INDL-1,
     $                  W(WBEGIN), WGAP(WBEGIN), WERR(WBEGIN),
     $                  WORK( 2*N+1 ), IWORK, PIVMIN, SPDIAM,
     $                  IN, IINFO )
            IF( IINFO .NE. 0 ) THEN
               INFO = -4
               RETURN
            END IF
<span class="comment">*</span><span class="comment">           <a name="SLARRB.668"></a><a href="slarrb.f.html#SLARRB.1">SLARRB</a> computes all gaps correctly except for the last one
</span><span class="comment">*</span><span class="comment">           Record distance to VU/GU
</span>            WGAP( WEND ) = MAX( ZERO,
     $           ( VU-SIGMA ) - ( W( WEND ) + WERR( WEND ) ) )
            DO 138 I = INDL, INDU
               M = M + 1
               IBLOCK(M) = JBLK
               INDEXW(M) = I
 138        CONTINUE
         ELSE
<span class="comment">*</span><span class="comment">           Call dqds to get all eigs (and then possibly delete unwanted
</span><span class="comment">*</span><span class="comment">           eigenvalues).
</span><span class="comment">*</span><span class="comment">           Note that dqds finds the eigenvalues of the L D L^T representation
</span><span class="comment">*</span><span class="comment">           of T to high relative accuracy. High relative accuracy
</span><span class="comment">*</span><span class="comment">           might be lost when the shift of the RRR is subtracted to obtain
</span><span class="comment">*</span><span class="comment">           the eigenvalues of T. However, T is not guaranteed to define its
</span><span class="comment">*</span><span class="comment">           eigenvalues to high relative accuracy anyway.
</span><span class="comment">*</span><span class="comment">           Set RTOL to the order of the tolerance used in <a name="SLASQ2.685"></a><a href="slasq2.f.html#SLASQ2.1">SLASQ2</a>
</span><span class="comment">*</span><span class="comment">           This is an ESTIMATED error, the worst case bound is 4*N*EPS
</span><span class="comment">*</span><span class="comment">           which is usually too large and requires unnecessary work to be
</span><span class="comment">*</span><span class="comment">           done by bisection when computing the eigenvectors
</span>            RTOL = LOG(REAL(IN)) * FOUR * EPS
            J = IBEGIN
            DO 140 I = 1, IN - 1
               WORK( 2*I-1 ) = ABS( D( J ) )
               WORK( 2*I ) = E( J )*E( J )*WORK( 2*I-1 )
               J = J + 1
  140       CONTINUE
            WORK( 2*IN-1 ) = ABS( D( IEND ) )
            WORK( 2*IN ) = ZERO
            CALL <a name="SLASQ2.698"></a><a href="slasq2.f.html#SLASQ2.1">SLASQ2</a>( IN, WORK, IINFO )
            IF( IINFO .NE. 0 ) THEN
<span class="comment">*</span><span class="comment">              If IINFO = -5 then an index is part of a tight cluster
</span><span class="comment">*</span><span class="comment">              and should be changed. The index is in IWORK(1) and the
</span><span class="comment">*</span><span class="comment">              gap is in WORK(N+1)
</span>               INFO = -5
               RETURN
            ELSE
<span class="comment">*</span><span class="comment">              Test that all eigenvalues are positive as expected
</span>               DO 149 I = 1, IN
                  IF( WORK( I ).LT.ZERO ) THEN
                     INFO = -6
                     RETURN
                  ENDIF
 149           CONTINUE
            END IF
            IF( SGNDEF.GT.ZERO ) THEN
               DO 150 I = INDL, INDU
                  M = M + 1
                  W( M ) = WORK( IN-I+1 )
                  IBLOCK( M ) = JBLK
                  INDEXW( M ) = I
 150           CONTINUE
            ELSE
               DO 160 I = INDL, INDU
                  M = M + 1
                  W( M ) = -WORK( I )
                  IBLOCK( M ) = JBLK
                  INDEXW( M ) = I
 160           CONTINUE
            END IF

            DO 165 I = M - MB + 1, M
<span class="comment">*</span><span class="comment">              the value of RTOL below should be the tolerance in <a name="SLASQ2.731"></a><a href="slasq2.f.html#SLASQ2.1">SLASQ2</a>
</span>               WERR( I ) = RTOL * ABS( W(I) )
 165        CONTINUE
            DO 166 I = M - MB + 1, M - 1
<span class="comment">*</span><span class="comment">              compute the right gap between the intervals
</span>               WGAP( I ) = MAX( ZERO,
     $                          W(I+1)-WERR(I+1) - (W(I)+WERR(I)) )
 166        CONTINUE
            WGAP( M ) = MAX( ZERO,
     $           ( VU-SIGMA ) - ( W( M ) + WERR( M ) ) )
         END IF
<span class="comment">*</span><span class="comment">        proceed with next block
</span>         IBEGIN = IEND + 1
         WBEGIN = WEND + 1
 170  CONTINUE
<span class="comment">*</span><span class="comment">
</span>
      RETURN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     end of <a name="SLARRE.750"></a><a href="slarre.f.html#SLARRE.1">SLARRE</a>
</span><span class="comment">*</span><span class="comment">
</span>      END

</pre>

 </body>
</html>
